module CanopyFluxesMod

!------------------------------------------------------------------------------
!BOP
!
! !MODULE: CanopyFluxesMod
!
! !DESCRIPTION:
! Calculates the leaf temperature and the leaf fluxes,
! transpiration, photosynthesis and  updates the dew
! accumulation due to evaporation.
!
! !USES:
   use abortutils,   only: endrun
   use clm_varctl,   only: iulog
   use shr_sys_mod,  only: shr_sys_flush
!
! !PUBLIC TYPES:
   implicit none
   save
!
! !PUBLIC MEMBER FUNCTIONS:
   public :: CanopyFluxes !Calculates the leaf temperature and leaf fluxes
!
! !PRIVATE MEMBER FUNCTIONS:
   private :: Stomata     !Leaf stomatal resistance and leaf photosynthesis
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
! 4/25/05, Peter Thornton: replaced old Stomata subroutine with what
!   used to be called StomataCN, as part of migration to new sun/shade
!   algorithms. 
!
!EOP
!------------------------------------------------------------------------------

contains

!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: CanopyFluxes 
!
! !INTERFACE:
  subroutine CanopyFluxes(lbg, ubg, lbc, ubc, lbp, ubp, &
                          num_nolakep, filter_nolakep)
!
! !DESCRIPTION:
! 1. Calculates the leaf temperature:
! 2. Calculates the leaf fluxes, transpiration, photosynthesis and
!    updates the dew accumulation due to evaporation.
!
! Method:
! Use the Newton-Raphson iteration to solve for the foliage
! temperature that balances the surface energy budget:
!
! f(t_veg) = Net radiation - Sensible - Latent = 0
! f(t_veg) + d(f)/d(t_veg) * dt_veg = 0     (*)
!
! Note:
! (1) In solving for t_veg, t_grnd is given from the previous timestep.
! (2) The partial derivatives of aerodynamical resistances, which cannot
!     be determined analytically, are ignored for d(H)/dT and d(LE)/dT
! (3) The weighted stomatal resistance of sunlit and shaded foliage is used
! (4) Canopy air temperature and humidity are derived from => Hc + Hg = Ha
!                                                          => Ec + Eg = Ea
! (5) Energy loss is due to: numerical truncation of energy budget equation
!     (*); and "ecidif" (see the code) which is dropped into the sensible
!     heat
! (6) The convergence criteria: the difference, del = t_veg(n+1)-t_veg(n)
!     and del2 = t_veg(n)-t_veg(n-1) less than 0.01 K, and the difference
!     of water flux from the leaf between the iteration step (n+1) and (n)
!     less than 0.1 W/m2; or the iterative steps over 40.
!
! !USES:
    use shr_kind_mod       , only : r8 => shr_kind_r8
    use clmtype
    use clm_atmlnd         , only : clm_a2l
    use clm_time_manager   , only : get_step_size, get_prev_date
    use clm_varpar         , only : nlevgrnd, nlevsno
    use clm_varcon         , only : sb, cpair, hvap, vkc, grav, denice, &
                                    denh2o, tfrz, csoilc, tlsai_crit, alpha_aero, &
                                    isecspday, degpsec
    use shr_const_mod      , only : SHR_CONST_TKFRZ
    use pftvarcon          , only : nirrig
    use QSatMod            , only : QSat
    use FrictionVelocityMod, only : FrictionVelocity, MoninObukIni
    use spmdMod            , only : masterproc
!
! !ARGUMENTS:
    implicit none
    integer, intent(in) :: lbg, ubg                    ! gridcell bounds
    integer, intent(in) :: lbc, ubc                    ! column bounds
    integer, intent(in) :: lbp, ubp                    ! pft bounds
    integer, intent(in) :: num_nolakep                 ! number of column non-lake points in pft filter
    integer, intent(in) :: filter_nolakep(ubp-lbp+1)   ! pft filter for non-lake points
!
! !CALLED FROM:
! subroutine Biogeophysics1 in module Biogeophysics1Mod
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
! 12/19/01, Peter Thornton
! Changed tg to t_grnd for consistency with other routines
! 1/29/02, Peter Thornton
! Migrate to new data structures, new calling protocol. For now co2 and
! o2 partial pressures are hardwired, but they should be coming in from
! forc_pco2 and forc_po2. Keeping the same hardwired values as in CLM2 to
! assure bit-for-bit results in the first comparisons.
! 27 February 2008: Keith Oleson; Sparse/dense aerodynamic parameters from
! X. Zeng
! 6 March 2009: Peter Thornton; Daylength control on Vcmax, from Bill Bauerle
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in variables
!
   integer , pointer :: frac_veg_nosno(:) ! frac of veg not covered by snow (0 OR 1 now) [-]
   integer , pointer :: ivt(:)         ! pft vegetation type
   integer , pointer :: pcolumn(:)     ! pft's column index
   integer , pointer :: plandunit(:)   ! pft's landunit index
   integer , pointer :: pgridcell(:)   ! pft's gridcell index
   real(r8), pointer :: forc_th(:)     ! atmospheric potential temperature (Kelvin)
   real(r8), pointer :: t_grnd(:)      ! ground surface temperature [K]
   real(r8), pointer :: thm(:)         ! intermediate variable (forc_t+0.0098*forc_hgt_t_pft)
   real(r8), pointer :: qg(:)          ! specific humidity at ground surface [kg/kg]
   real(r8), pointer :: thv(:)         ! virtual potential temperature (kelvin)
   real(r8), pointer :: z0mv(:)        ! roughness length over vegetation, momentum [m]
   real(r8), pointer :: z0hv(:)        ! roughness length over vegetation, sensible heat [m]
   real(r8), pointer :: z0qv(:)        ! roughness length over vegetation, latent heat [m]
   real(r8), pointer :: z0mg(:)        ! roughness length of ground, momentum [m]
   real(r8), pointer :: dqgdT(:)       ! temperature derivative of "qg"
   real(r8), pointer :: htvp(:)        ! latent heat of evaporation (/sublimation) [J/kg]
   real(r8), pointer :: emv(:)         ! ground emissivity
   real(r8), pointer :: emg(:)         ! vegetation emissivity
   real(r8), pointer :: forc_pbot(:)   ! atmospheric pressure (Pa)
   real(r8), pointer :: forc_pco2(:)   ! partial pressure co2 (Pa)
#if (defined C13)
   ! 4/14/05: PET
   ! Adding isotope code
   real(r8), pointer :: forc_pc13o2(:) ! partial pressure c13o2 (Pa)
#endif
   
   real(r8), pointer :: forc_po2(:)    ! partial pressure o2 (Pa)
   real(r8), pointer :: forc_q(:)      ! atmospheric specific humidity (kg/kg)
   real(r8), pointer :: forc_u(:)      ! atmospheric wind speed in east direction (m/s)
   real(r8), pointer :: forc_v(:)      ! atmospheric wind speed in north direction (m/s)
   real(r8), pointer :: forc_hgt_u_pft(:) !observational height of wind at pft level [m]
   real(r8), pointer :: forc_rho(:)    ! density (kg/m**3)
   real(r8), pointer :: forc_lwrad(:)  ! downward infrared (longwave) radiation (W/m**2)
   real(r8), pointer :: displa(:)      ! displacement height (m)
   real(r8), pointer :: elai(:)        ! one-sided leaf area index with burying by snow
   real(r8), pointer :: esai(:)        ! one-sided stem area index with burying by snow
   real(r8), pointer :: fdry(:)        ! fraction of foliage that is green and dry [-]
   real(r8), pointer :: fwet(:)        ! fraction of canopy that is wet (0 to 1)
   real(r8), pointer :: laisun(:)      ! sunlit leaf area
   real(r8), pointer :: laisha(:)      ! shaded leaf area
   real(r8), pointer :: sabv(:)        ! solar radiation absorbed by vegetation (W/m**2)
   real(r8), pointer :: watsat(:,:)    ! volumetric soil water at saturation (porosity)
   real(r8), pointer :: watdry(:,:)    ! btran parameter for btran=0
   real(r8), pointer :: watopt(:,:)    ! btran parameter for btran = 1
   real(r8), pointer :: h2osoi_ice(:,:)! ice lens (kg/m2)
   real(r8), pointer :: h2osoi_liq(:,:)! liquid water (kg/m2)
   real(r8), pointer :: dz(:,:)        ! layer depth (m)
   real(r8), pointer :: t_soisno(:,:)  ! soil temperature (Kelvin)
   real(r8), pointer :: sucsat(:,:)    ! minimum soil suction (mm)
   real(r8), pointer :: bsw(:,:)       ! Clapp and Hornberger "b"
   real(r8), pointer :: rootfr(:,:)    ! fraction of roots in each soil layer
   real(r8), pointer :: dleaf(:)       ! characteristic leaf dimension (m)
   real(r8), pointer :: smpso(:)       ! soil water potential at full stomatal opening (mm)
   real(r8), pointer :: smpsc(:)       ! soil water potential at full stomatal closure (mm)
   real(r8), pointer :: frac_sno(:)    ! fraction of ground covered by snow (0 to 1)
   real(r8), pointer :: htop(:)        ! canopy top(m)
   real(r8), pointer :: snowdp(:)      ! snow height (m)
   real(r8), pointer :: soilbeta(:)    ! soil wetness relative to field capacity
   real(r8), pointer :: lat(:)         ! latitude (radians)
   real(r8), pointer :: decl(:)        ! declination angle (radians)
   real(r8), pointer :: max_dayl(:)    !maximum daylength for this column (s)
   real(r8), pointer :: londeg(:)      ! longitude (degrees) (for calculation of local time)
   
!
! local pointers to implicit inout arguments
!
   real(r8), pointer :: cgrnds(:)      ! deriv. of soil sensible heat flux wrt soil temp [w/m2/k]
   real(r8), pointer :: cgrndl(:)      ! deriv. of soil latent heat flux wrt soil temp [w/m**2/k]
   real(r8), pointer :: t_veg(:)       ! vegetation temperature (Kelvin)
   real(r8), pointer :: t_ref2m(:)     ! 2 m height surface air temperature (Kelvin)
   real(r8), pointer :: q_ref2m(:)     ! 2 m height surface specific humidity (kg/kg)
   real(r8), pointer :: t_ref2m_r(:)   ! Rural 2 m height surface air temperature (Kelvin)
   real(r8), pointer :: rh_ref2m(:)    ! 2 m height surface relative humidity (%)
   real(r8), pointer :: rh_ref2m_r(:)  ! Rural 2 m height surface relative humidity (%)
   real(r8), pointer :: h2ocan(:)      ! canopy water (mm H2O)
   real(r8), pointer :: cisun(:)       !sunlit intracellular CO2 (Pa)
   real(r8), pointer :: cisha(:)       !shaded intracellular CO2 (Pa)
!
! local pointers to implicit out arguments
!
   real(r8), pointer :: rb1(:)             ! boundary layer resistance (s/m)
   real(r8), pointer :: cgrnd(:)           ! deriv. of soil energy flux wrt to soil temp [w/m2/k]
   real(r8), pointer :: dlrad(:)           ! downward longwave radiation below the canopy [W/m2]
   real(r8), pointer :: ulrad(:)           ! upward longwave radiation above the canopy [W/m2]
   real(r8), pointer :: ram1(:)            ! aerodynamical resistance (s/m)
   real(r8), pointer :: btran(:)           ! transpiration wetness factor (0 to 1)
   real(r8), pointer :: rssun(:)           ! sunlit stomatal resistance (s/m)
   real(r8), pointer :: rssha(:)           ! shaded stomatal resistance (s/m)
   real(r8), pointer :: psnsun(:)          ! sunlit leaf photosynthesis (umol CO2 /m**2/ s)
   real(r8), pointer :: psnsha(:)          ! shaded leaf photosynthesis (umol CO2 /m**2/ s)
#if (defined C13)
   ! 4/14/05: PET
   ! Adding isotope code
   real(r8), pointer :: c13_psnsun(:)      ! sunlit leaf photosynthesis (umol 13CO2 /m**2/ s)
   real(r8), pointer :: c13_psnsha(:)      ! shaded leaf photosynthesis (umol 13CO2 /m**2/ s)
   ! 4/21/05: PET
   ! Adding isotope code
	real(r8), pointer :: rc13_canair(:)     !C13O2/C12O2 in canopy air
   real(r8), pointer :: rc13_psnsun(:)     !C13O2/C12O2 in sunlit canopy psn flux
   real(r8), pointer :: rc13_psnsha(:)     !C13O2/C12O2 in shaded canopy psn flux
   real(r8), pointer :: alphapsnsun(:)     !fractionation factor in sunlit canopy psn flux
   real(r8), pointer :: alphapsnsha(:)     !fractionation factor in shaded canopy psn flux
#endif
   
   real(r8), pointer :: qflx_tran_veg(:)   ! vegetation transpiration (mm H2O/s) (+ = to atm)
   real(r8), pointer :: dt_veg(:)          ! change in t_veg, last iteration (Kelvin)
   real(r8), pointer :: qflx_evap_veg(:)   ! vegetation evaporation (mm H2O/s) (+ = to atm)
   real(r8), pointer :: eflx_sh_veg(:)     ! sensible heat flux from leaves (W/m**2) [+ to atm]
   real(r8), pointer :: taux(:)            ! wind (shear) stress: e-w (kg/m/s**2)
   real(r8), pointer :: tauy(:)            ! wind (shear) stress: n-s (kg/m/s**2)
   real(r8), pointer :: eflx_sh_grnd(:)    ! sensible heat flux from ground (W/m**2) [+ to atm]
   real(r8), pointer :: qflx_evap_soi(:)   ! soil evaporation (mm H2O/s) (+ = to atm)
   real(r8), pointer :: fpsn(:)            ! photosynthesis (umol CO2 /m**2 /s)
   real(r8), pointer :: rootr(:,:)         ! effective fraction of roots in each soil layer
   real(r8), pointer :: rresis(:,:)        ! root resistance by layer (0-1)  (nlevgrnd)	
   real(r8), pointer :: irrig_rate(:)      ! current irrigation rate [mm/s]
   integer, pointer  :: n_irrig_steps_left(:) ! # of time steps for which we still need to irrigate today
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
   real(r8), parameter :: btran0 = 0.0_r8  ! initial value
   real(r8), parameter :: zii = 1000.0_r8  ! convective boundary layer height [m]
   real(r8), parameter :: beta = 1.0_r8    ! coefficient of conective velocity [-]
   real(r8), parameter :: delmax = 1.0_r8  ! maxchange in  leaf temperature [K]
   real(r8), parameter :: dlemin = 0.1_r8  ! max limit for energy flux convergence [w/m2]
   real(r8), parameter :: dtmin = 0.01_r8  ! max limit for temperature convergence [K]
   integer , parameter :: itmax = 40       ! maximum number of iteration [-]
   integer , parameter :: itmin = 2        ! minimum number of iteration [-]
   real(r8), parameter :: irrig_min_lai = 0.0_r8           ! Minimum LAI for irrigation
   real(r8), parameter :: irrig_btran_thresh = 0.999999_r8 ! Irrigate when btran falls below 0.999999 rather than 1 to allow for round-off error
   integer , parameter :: irrig_start_time = isecspday/4   ! Time of day to check whether we need irrigation, seconds (0 = midnight). 
                                                           ! We start applying the irrigation in the time step FOLLOWING this time, 
                                                           ! since we won't begin irrigating until the next call to Hydrology1
   integer , parameter :: irrig_length = isecspday/6       ! Desired amount of time to irrigate per day (sec). Actual time may 
                                                           ! differ if this is not a multiple of dtime. Irrigation won't work properly 
                                                           ! if dtime > secsperday
   real(r8), parameter :: irrig_factor = 0.7_r8            ! Determines target soil moisture level for irrigation. If h2osoi_liq_so 
                                                           ! is the soil moisture level at which stomata are fully open and 
                                                           ! h2osoi_liq_sat is the soil moisture level at saturation (eff_porosity), 
                                                           ! then the target soil moisture level is 
                                                           !     (h2osoi_liq_so + irrig_factor*(h2osoi_liq_sat - h2osoi_liq_so)). 
                                                           ! A value of 0 means that the target soil moisture level is h2osoi_liq_so. 
                                                           ! A value of 1 means that the target soil moisture level is h2osoi_liq_sat

   !added by K.Sakaguchi for litter resistance
   real(r8), parameter :: lai_dl = 0.5_r8  ! placeholder for (dry) plant litter area index (m2/m2)
   real(r8), parameter :: z_dl = 0.05_r8   ! placeholder for (dry) litter layer thickness (m)
   !added by K.Sakaguchi for stability formulation
   real(r8), parameter :: ria  = 0.5_r8    ! free parameter for stable formulation (currently = 0.5, "gamma" in Sakaguchi&Zeng,2008)
   real(r8) :: dtime                 ! land model time step (sec)
   real(r8) :: zldis(lbp:ubp)        ! reference height "minus" zero displacement height [m]
   real(r8) :: zeta                  ! dimensionless height used in Monin-Obukhov theory
   real(r8) :: wc                    ! convective velocity [m/s]
   real(r8) :: dth(lbp:ubp)          ! diff of virtual temp. between ref. height and surface
   real(r8) :: dthv(lbp:ubp)         ! diff of vir. poten. temp. between ref. height and surface
   real(r8) :: dqh(lbp:ubp)          ! diff of humidity between ref. height and surface
   real(r8) :: obu(lbp:ubp)          ! Monin-Obukhov length (m)
   real(r8) :: um(lbp:ubp)           ! wind speed including the stablity effect [m/s]
   real(r8) :: ur(lbp:ubp)           ! wind speed at reference height [m/s]
   real(r8) :: uaf(lbp:ubp)          ! velocity of air within foliage [m/s]
   real(r8) :: temp1(lbp:ubp)        ! relation for potential temperature profile
   real(r8) :: temp12m(lbp:ubp)      ! relation for potential temperature profile applied at 2-m
   real(r8) :: temp2(lbp:ubp)        ! relation for specific humidity profile
   real(r8) :: temp22m(lbp:ubp)      ! relation for specific humidity profile applied at 2-m
   real(r8) :: ustar(lbp:ubp)        ! friction velocity [m/s]
   real(r8) :: tstar                 ! temperature scaling parameter
   real(r8) :: qstar                 ! moisture scaling parameter
   real(r8) :: thvstar               ! virtual potential temperature scaling parameter
   real(r8) :: taf(lbp:ubp)          ! air temperature within canopy space [K]
   real(r8) :: qaf(lbp:ubp)          ! humidity of canopy air [kg/kg]
   real(r8) :: rpp                   ! fraction of potential evaporation from leaf [-]
   real(r8) :: rppdry                ! fraction of potential evaporation through transp [-]
   real(r8) :: cf                    ! heat transfer coefficient from leaves [-]
   real(r8) :: rb(lbp:ubp)           ! leaf boundary layer resistance [s/m]
   real(r8) :: rah(lbp:ubp,2)        ! thermal resistance [s/m]
   real(r8) :: raw(lbp:ubp,2)        ! moisture resistance [s/m]
   real(r8) :: wta                   ! heat conductance for air [m/s]
   real(r8) :: wtg(lbp:ubp)          ! heat conductance for ground [m/s]
   real(r8) :: wtl                   ! heat conductance for leaf [m/s]
   real(r8) :: wta0(lbp:ubp)         ! normalized heat conductance for air [-]
   real(r8) :: wtl0(lbp:ubp)         ! normalized heat conductance for leaf [-]
   real(r8) :: wtg0                  ! normalized heat conductance for ground [-]
   real(r8) :: wtal(lbp:ubp)         ! normalized heat conductance for air and leaf [-]
   real(r8) :: wtga                  ! normalized heat cond. for air and ground  [-]
   real(r8) :: wtaq                  ! latent heat conductance for air [m/s]
   real(r8) :: wtlq                  ! latent heat conductance for leaf [m/s]
   real(r8) :: wtgq(lbp:ubp)         ! latent heat conductance for ground [m/s]
   real(r8) :: wtaq0(lbp:ubp)        ! normalized latent heat conductance for air [-]
   real(r8) :: wtlq0(lbp:ubp)        ! normalized latent heat conductance for leaf [-]
   real(r8) :: wtgq0                 ! normalized heat conductance for ground [-]
   real(r8) :: wtalq(lbp:ubp)        ! normalized latent heat cond. for air and leaf [-]
   real(r8) :: wtgaq                 ! normalized latent heat cond. for air and ground [-]
   real(r8) :: el(lbp:ubp)           ! vapor pressure on leaf surface [pa]
   real(r8) :: deldT                 ! derivative of "el" on "t_veg" [pa/K]
   real(r8) :: qsatl(lbp:ubp)        ! leaf specific humidity [kg/kg]
   real(r8) :: qsatldT(lbp:ubp)      ! derivative of "qsatl" on "t_veg"
   real(r8) :: e_ref2m               ! 2 m height surface saturated vapor pressure [Pa]
   real(r8) :: de2mdT                ! derivative of 2 m height surface saturated vapor pressure on t_ref2m
   real(r8) :: qsat_ref2m            ! 2 m height surface saturated specific humidity [kg/kg]
   real(r8) :: dqsat2mdT             ! derivative of 2 m height surface saturated specific humidity on t_ref2m
   real(r8) :: air(lbp:ubp),bir(lbp:ubp),cir(lbp:ubp)  ! atmos. radiation temporay set
   real(r8) :: dc1,dc2               ! derivative of energy flux [W/m2/K]
   real(r8) :: delt                  ! temporary
   real(r8) :: delq(lbp:ubp)         ! temporary
   real(r8) :: del(lbp:ubp)          ! absolute change in leaf temp in current iteration [K]
   real(r8) :: del2(lbp:ubp)         ! change in leaf temperature in previous iteration [K]
   real(r8) :: dele(lbp:ubp)         ! change in latent heat flux from leaf [K]
   real(r8) :: dels                  ! change in leaf temperature in current iteration [K]
   real(r8) :: det(lbp:ubp)          ! maximum leaf temp. change in two consecutive iter [K]
   real(r8) :: efeb(lbp:ubp)         ! latent heat flux from leaf (previous iter) [mm/s]
   real(r8) :: efeold                ! latent heat flux from leaf (previous iter) [mm/s]
   real(r8) :: efpot                 ! potential latent energy flux [kg/m2/s]
   real(r8) :: efe(lbp:ubp)          ! water flux from leaf [mm/s]
   real(r8) :: efsh                  ! sensible heat from leaf [mm/s]
   real(r8) :: obuold(lbp:ubp)       ! monin-obukhov length from previous iteration
   real(r8) :: tlbef(lbp:ubp)        ! leaf temperature from previous iteration [K]
   real(r8) :: ecidif                ! excess energies [W/m2]
   real(r8) :: err(lbp:ubp)          ! balance error
   real(r8) :: erre                  ! balance error
   real(r8) :: co2(lbp:ubp)          ! atmospheric co2 partial pressure (pa)
#if (defined C13)
   ! 4/14/05: PET
   ! Adding isotope code
   real(r8) :: c13o2(lbp:ubp)        ! atmospheric c13o2 partial pressure (pa)
#endif
   
   real(r8) :: o2(lbp:ubp)           ! atmospheric o2 partial pressure (pa)
   real(r8) :: svpts(lbp:ubp)        ! saturation vapor pressure at t_veg (pa)
   real(r8) :: eah(lbp:ubp)          ! canopy air vapor pressure (pa)
   real(r8) :: s_node                ! vol_liq/eff_porosity
   real(r8) :: smp_node              ! matrix potential
   real(r8) :: vol_ice               ! partial volume of ice lens in layer
   real(r8) :: eff_porosity          ! effective porosity in layer
   real(r8) :: vol_liq               ! partial volume of liquid water in layer
   integer  :: itlef                 ! counter for leaf temperature iteration [-]
   integer  :: nmozsgn(lbp:ubp)      ! number of times stability changes sign
   real(r8) :: w                     ! exp(-LSAI)
   real(r8) :: csoilcn               ! interpolated csoilc for less than dense canopies
   real(r8) :: fm(lbp:ubp)           ! needed for BGC only to diagnose 10m wind speed
   real(r8) :: wtshi                 ! sensible heat resistance for air, grnd and leaf [-]
   real(r8) :: wtsqi                 ! latent heat resistance for air, grnd and leaf [-]
   integer  :: j                     ! soil/snow level index
   integer  :: p                     ! pft index
   integer  :: c                     ! column index
   integer  :: l                     ! landunit index
   integer  :: g                     ! gridcell index
   integer  :: fp                    ! lake filter pft index
   integer  :: fn                    ! number of values in pft filter
   integer  :: fnorig                ! number of values in pft filter copy
   integer  :: fnold                 ! temporary copy of pft count
   integer  :: f                     ! filter index
   integer  :: filterp(ubp-lbp+1)    ! temporary filter
   integer  :: fporig(ubp-lbp+1)     ! temporary filter
   real(r8) :: displa_loc(lbp:ubp)   ! temporary copy
   real(r8) :: z0mv_loc(lbp:ubp)     ! temporary copy
   real(r8) :: z0hv_loc(lbp:ubp)     ! temporary copy
   real(r8) :: z0qv_loc(lbp:ubp)     ! temporary copy
   logical  :: found                 ! error flag for canopy above forcing hgt
   integer  :: index                 ! pft index for error
   real(r8) :: egvf                  ! effective green vegetation fraction
   real(r8) :: lt                    ! elai+esai
   real(r8) :: ri                    ! stability parameter for under canopy air (unitless)
   real(r8) :: csoilb                ! turbulent transfer coefficient over bare soil (unitless)
   real(r8) :: ricsoilc              ! modified transfer coefficient under dense canopy (unitless)
   real(r8) :: snowdp_c              ! critical snow depth to cover plant litter (m)
   real(r8) :: rdl                   ! dry litter layer resistance for water vapor  (s/m)
   real(r8) :: elai_dl               ! exposed (dry) plant litter area index
   real(r8) :: fsno_dl               ! effective snow cover over plant litter
   real(r8) :: dayl                  ! daylength (s)
   real(r8) :: temp                  ! temporary, for daylength calculation
   real(r8) :: dayl_factor(lbp:ubp)  ! scalar (0-1) for daylength effect on Vcmax
   integer  :: yr                       ! year at start of time step
   integer  :: mon                      ! month at start of time step
   integer  :: day                      ! day at start of time step
   integer  :: time                     ! time at start of time step (seconds after 0Z)
   integer  :: local_time               ! local time at start of time step (seconds after solar midnight)
   integer  :: irrig_nsteps_per_day     ! number of time steps per day in which we irrigate
   logical  :: check_for_irrig(lbp:ubp) ! where do we need to check soil moisture to see if we need to irrigate?
   logical  :: frozen_soil(lbc:ubc)     ! set to true if we have encountered a frozen soil layer
   real(r8) :: vol_liq_so               ! partial volume of liquid water in layer for which smp_node = smpso
   real(r8) :: h2osoi_liq_so            ! liquid water corresponding to vol_liq_so for this layer [kg/m2]
   real(r8) :: h2osoi_liq_sat           ! liquid water corresponding to eff_porosity for this layer [kg/m2]
   real(r8) :: deficit                  ! difference between desired soil moisture level for this layer and current soil moisture level [kg/m2]
!------------------------------------------------------------------------------

   ! Assign local pointers to derived type members (gridcell-level)

   forc_lwrad     => clm_a2l%forc_lwrad
   forc_pco2      => clm_a2l%forc_pco2
#if (defined C13)
   forc_pc13o2    => clm_a2l%forc_pc13o2
#endif
   forc_po2       => clm_a2l%forc_po2
   forc_q         => clm_a2l%forc_q
   forc_pbot      => clm_a2l%forc_pbot
   forc_u         => clm_a2l%forc_u
   forc_v         => clm_a2l%forc_v
   forc_th        => clm_a2l%forc_th
   forc_rho       => clm_a2l%forc_rho
   lat            => clm3%g%lat
   londeg         => clm3%g%londeg

   ! Assign local pointers to derived type members (column-level)

   t_soisno       => clm3%g%l%c%ces%t_soisno
   watsat         => clm3%g%l%c%cps%watsat
   watdry         => clm3%g%l%c%cps%watdry 
   watopt         => clm3%g%l%c%cps%watopt 
   h2osoi_ice     => clm3%g%l%c%cws%h2osoi_ice
   dz             => clm3%g%l%c%cps%dz
   h2osoi_liq     => clm3%g%l%c%cws%h2osoi_liq
   sucsat         => clm3%g%l%c%cps%sucsat
   bsw            => clm3%g%l%c%cps%bsw
   emg            => clm3%g%l%c%cps%emg
   t_grnd         => clm3%g%l%c%ces%t_grnd
   qg             => clm3%g%l%c%cws%qg
   thv            => clm3%g%l%c%ces%thv
   dqgdT          => clm3%g%l%c%cws%dqgdT
   htvp           => clm3%g%l%c%cps%htvp
   z0mg           => clm3%g%l%c%cps%z0mg
   frac_sno       => clm3%g%l%c%cps%frac_sno
   snowdp         => clm3%g%l%c%cps%snowdp
   soilbeta       => clm3%g%l%c%cws%soilbeta
   decl           => clm3%g%l%c%cps%decl
   max_dayl       => clm3%g%l%c%cps%max_dayl

   ! Assign local pointers to derived type members (pft-level)

   rb1            => clm3%g%l%c%p%pps%rb1
   ivt            => clm3%g%l%c%p%itype
   pcolumn        => clm3%g%l%c%p%column
   plandunit      => clm3%g%l%c%p%landunit
   pgridcell      => clm3%g%l%c%p%gridcell
   frac_veg_nosno => clm3%g%l%c%p%pps%frac_veg_nosno
   btran          => clm3%g%l%c%p%pps%btran
   rootfr         => clm3%g%l%c%p%pps%rootfr
   rootr          => clm3%g%l%c%p%pps%rootr
   rresis         => clm3%g%l%c%p%pps%rresis
   emv            => clm3%g%l%c%p%pps%emv
   t_veg          => clm3%g%l%c%p%pes%t_veg
   displa         => clm3%g%l%c%p%pps%displa
   z0mv           => clm3%g%l%c%p%pps%z0mv
   z0hv           => clm3%g%l%c%p%pps%z0hv
   z0qv           => clm3%g%l%c%p%pps%z0qv
   ram1           => clm3%g%l%c%p%pps%ram1
   htop           => clm3%g%l%c%p%pps%htop
   rssun          => clm3%g%l%c%p%pps%rssun
   rssha          => clm3%g%l%c%p%pps%rssha
   cisun          => clm3%g%l%c%p%pps%cisun
   cisha          => clm3%g%l%c%p%pps%cisha
   psnsun         => clm3%g%l%c%p%pcf%psnsun
   psnsha         => clm3%g%l%c%p%pcf%psnsha
#if (defined C13)
   ! 4/14/05: PET
   ! Adding isotope code
   c13_psnsun     => clm3%g%l%c%p%pc13f%psnsun
   c13_psnsha     => clm3%g%l%c%p%pc13f%psnsha
   ! 4/21/05: PET
   ! Adding isotope code
   rc13_canair    => clm3%g%l%c%p%pepv%rc13_canair
   rc13_psnsun    => clm3%g%l%c%p%pepv%rc13_psnsun
   rc13_psnsha    => clm3%g%l%c%p%pepv%rc13_psnsha
   alphapsnsun    => clm3%g%l%c%p%pps%alphapsnsun
   alphapsnsha    => clm3%g%l%c%p%pps%alphapsnsha
#endif
   
   elai           => clm3%g%l%c%p%pps%elai
   esai           => clm3%g%l%c%p%pps%esai
   fdry           => clm3%g%l%c%p%pps%fdry
   laisun         => clm3%g%l%c%p%pps%laisun
   laisha         => clm3%g%l%c%p%pps%laisha
   qflx_tran_veg  => clm3%g%l%c%p%pwf%qflx_tran_veg
   fwet           => clm3%g%l%c%p%pps%fwet
   h2ocan         => clm3%g%l%c%p%pws%h2ocan
   dt_veg         => clm3%g%l%c%p%pps%dt_veg
   sabv           => clm3%g%l%c%p%pef%sabv
   qflx_evap_veg  => clm3%g%l%c%p%pwf%qflx_evap_veg
   eflx_sh_veg    => clm3%g%l%c%p%pef%eflx_sh_veg
   taux           => clm3%g%l%c%p%pmf%taux
   tauy           => clm3%g%l%c%p%pmf%tauy
   eflx_sh_grnd   => clm3%g%l%c%p%pef%eflx_sh_grnd
   qflx_evap_soi  => clm3%g%l%c%p%pwf%qflx_evap_soi
   t_ref2m        => clm3%g%l%c%p%pes%t_ref2m
   q_ref2m        => clm3%g%l%c%p%pes%q_ref2m
   t_ref2m_r      => clm3%g%l%c%p%pes%t_ref2m_r
   rh_ref2m_r     => clm3%g%l%c%p%pes%rh_ref2m_r
   rh_ref2m       => clm3%g%l%c%p%pes%rh_ref2m
   dlrad          => clm3%g%l%c%p%pef%dlrad
   ulrad          => clm3%g%l%c%p%pef%ulrad
   cgrnds         => clm3%g%l%c%p%pef%cgrnds
   cgrndl         => clm3%g%l%c%p%pef%cgrndl
   cgrnd          => clm3%g%l%c%p%pef%cgrnd
   fpsn           => clm3%g%l%c%p%pcf%fpsn
   forc_hgt_u_pft => clm3%g%l%c%p%pps%forc_hgt_u_pft
   thm            => clm3%g%l%c%p%pes%thm
   irrig_rate         => clm3%g%l%c%cps%irrig_rate
   n_irrig_steps_left => clm3%g%l%c%cps%n_irrig_steps_left
      
   ! Assign local pointers to derived type members (ecophysiological)

   dleaf          => pftcon%dleaf
   smpso          => pftcon%smpso
   smpsc          => pftcon%smpsc

   ! Determine step size

   dtime = get_step_size()
   irrig_nsteps_per_day = ((irrig_length + (dtime - 1))/dtime)  ! round up

   ! Filter pfts where frac_veg_nosno is non-zero

   fn = 0
   do fp = 1,num_nolakep
      p = filter_nolakep(fp)
      if (frac_veg_nosno(p) /= 0) then
         fn = fn + 1
         filterp(fn) = p
      end if
   end do

   ! Initialize

   do f = 1, fn
      p = filterp(f)
      del(p)    = 0._r8  ! change in leaf temperature from previous iteration
      efeb(p)   = 0._r8  ! latent head flux from leaf for previous iteration
      wtlq0(p)  = 0._r8
      wtalq(p)  = 0._r8
      wtgq(p)   = 0._r8
      wtaq0(p)  = 0._r8
      obuold(p) = 0._r8
      btran(p)  = btran0
   end do
   
   ! calculate daylength control for Vcmax
   do f = 1, fn
      p=filterp(f)
      c=pcolumn(p)
      g=pgridcell(p)
      ! calculate daylength
      temp = -(sin(lat(g))*sin(decl(c)))/(cos(lat(g)) * cos(decl(c)))
      temp = min(1._r8,max(-1._r8,temp))
      dayl = 2.0_r8 * 13750.9871_r8 * acos(temp)
      ! calculate dayl_factor as the ratio of (current:max dayl)^2
      ! set a minimum of 0.01 (1%) for the dayl_factor
      dayl_factor(p)=min(1._r8,max(0.01_r8,(dayl*dayl)/(max_dayl(c)*max_dayl(c))))
   end do

   rb1(lbp:ubp) = 0._r8

   ! Effective porosity of soil, partial volume of ice and liquid (needed for btran)
   ! and root resistance factors

   do j = 1,nlevgrnd
      do f = 1, fn
         p = filterp(f)
         c = pcolumn(p)
         l = plandunit(p)

         ! Root resistance factors

         vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice))
         eff_porosity = watsat(c,j)-vol_ice
         vol_liq = min(eff_porosity, h2osoi_liq(c,j)/(dz(c,j)*denh2o))
         if (vol_liq .le. 0._r8 .or. t_soisno(c,j) .le. tfrz-2._r8) then
            rootr(p,j) = 0._r8
         else
            s_node = max(vol_liq/eff_porosity,0.01_r8)
            smp_node = max(smpsc(ivt(p)), -sucsat(c,j)*s_node**(-bsw(c,j)))

            rresis(p,j) = min( (eff_porosity/watsat(c,j))* &
                          (smp_node - smpsc(ivt(p))) / (smpso(ivt(p)) - smpsc(ivt(p))), 1._r8)
            rootr(p,j) = rootfr(p,j)*rresis(p,j)
            btran(p) = btran(p) + rootr(p,j)
         endif 
      end do
   end do

   ! Normalize root resistances to get layer contribution to ET

   do j = 1,nlevgrnd
      do f = 1, fn
         p = filterp(f)
         if (btran(p) .gt. btran0) then
           rootr(p,j) = rootr(p,j)/btran(p)
         else
           rootr(p,j) = 0._r8
         end if
      end do
   end do

   ! Determine if irrigation is needed (over irrigated soil columns)

   ! First, determine in what grid cells we need to bother 'measuring' soil water, to see if we need irrigation
   ! Also set n_irrig_steps_left for these grid cells
   ! n_irrig_steps_left(p) > 0 is ok even if irrig_rate(p) ends up = 0
   ! in this case, we'll irrigate by 0 for the given number of time steps
   call get_prev_date(yr, mon, day, time)  ! get time as of beginning of time step
   do f = 1, fn
      p = filterp(f)
      c = pcolumn(p)
      g = pgridcell(p)
      if (ivt(p) == nirrig .and. elai(p) > irrig_min_lai .and. btran(p) < irrig_btran_thresh) then
         ! see if it's the right time of day to start irrigating:
         local_time = modulo(time + nint(londeg(g)/degpsec), isecspday)
         if (modulo(local_time - irrig_start_time, isecspday) < dtime) then
            ! it's time to start irrigating
            check_for_irrig(p)    = .true.
            n_irrig_steps_left(c) = irrig_nsteps_per_day
            irrig_rate(c)         = 0._r8  ! reset; we'll add to this later
         else
            check_for_irrig(p)    = .false.
         end if
      else  ! non-irrig pft or elai<=irrig_min_lai or btran>irrig_btran_thresh
         check_for_irrig(p)       = .false.
      end if

   end do


   ! Now 'measure' soil water for the grid cells identified above and see if the soil is dry enough to warrant irrigation
   frozen_soil(:) = .false.
   do j = 1,nlevgrnd
      do f = 1, fn
         p = filterp(f)
         c = pcolumn(p)
         if (check_for_irrig(p) .and. .not. frozen_soil(c)) then
            ! if level L was frozen, then we don't look at any levels below L
            if (t_soisno(c,j) <= SHR_CONST_TKFRZ) then
               frozen_soil(c) = .true.
            else if (rootfr(p,j) > 0._r8) then
               ! determine soil water deficit in this layer:

               ! Calculate vol_liq_so - i.e., vol_liq at which smp_node = smpso - by inverting the above equations 
               ! for the root resistance factors
               vol_ice      = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice))  ! this duplicates the above equation for vol_ice
               eff_porosity = watsat(c,j)-vol_ice  ! this duplicates the above equation for eff_porosity
               vol_liq_so   = eff_porosity * (-smpso(ivt(p))/sucsat(c,j))**(-1/bsw(c,j))

               ! Translate vol_liq_so and eff_porosity into h2osoi_liq_so and h2osoi_liq_sat and calculate deficit
               h2osoi_liq_so  = vol_liq_so * denh2o * dz(c,j)
               h2osoi_liq_sat = eff_porosity * denh2o * dz(c,j)
               deficit        = max((h2osoi_liq_so + irrig_factor*(h2osoi_liq_sat - h2osoi_liq_so)) - h2osoi_liq(c,j), 0._r8)

               ! Add deficit to irrig_rate, converting units from mm to mm/sec
               irrig_rate(c)  = irrig_rate(c) + deficit/(dtime*irrig_nsteps_per_day)

            end if  ! else if (rootfr(p,j) .gt. 0)
         end if     ! if (check_for_irrig(p) .and. .not. frozen_soil(c))
      end do        ! do f
   end do           ! do j

 ! Modify aerodynamic parameters for sparse/dense canopy (X. Zeng)
   do f = 1, fn
      p = filterp(f)
      c = pcolumn(p)

      lt = min(elai(p)+esai(p), tlsai_crit)
      egvf =(1._r8 - alpha_aero * exp(-lt)) / (1._r8 - alpha_aero * exp(-tlsai_crit))
      displa(p) = egvf * displa(p)
      z0mv(p)   = exp(egvf * log(z0mv(p)) + (1._r8 - egvf) * log(z0mg(c)))
      z0hv(p)   = z0mv(p)
      z0qv(p)   = z0mv(p)

  end do

   found = .false.
   do f = 1, fn
      p = filterp(f)
      c = pcolumn(p)
      g = pgridcell(p)

      ! Net absorbed longwave radiation by canopy and ground
      ! =air+bir*t_veg**4+cir*t_grnd(c)**4

      air(p) =   emv(p) * (1._r8+(1._r8-emv(p))*(1._r8-emg(c))) * forc_lwrad(g)
      bir(p) = - (2._r8-emv(p)*(1._r8-emg(c))) * emv(p) * sb
      cir(p) =   emv(p)*emg(c)*sb

      ! Saturated vapor pressure, specific humidity, and their derivatives
      ! at the leaf surface

      call QSat (t_veg(p), forc_pbot(g), el(p), deldT, qsatl(p), qsatldT(p))

      ! Determine atmospheric co2 and o2

      co2(p) = forc_pco2(g)
      o2(p)  = forc_po2(g)
      
#if (defined C13)
      ! 4/14/05: PET
      ! Adding isotope code
      c13o2(p) = forc_pc13o2(g)
#endif
      
      ! Initialize flux profile

      nmozsgn(p) = 0

      taf(p) = (t_grnd(c) + thm(p))/2._r8
      qaf(p) = (forc_q(g)+qg(c))/2._r8

      ur(p) = max(1.0_r8,sqrt(forc_u(g)*forc_u(g)+forc_v(g)*forc_v(g)))
      dth(p) = thm(p)-taf(p)
      dqh(p) = forc_q(g)-qaf(p)
      delq(p) = qg(c) - qaf(p)
      dthv(p) = dth(p)*(1._r8+0.61_r8*forc_q(g))+0.61_r8*forc_th(g)*dqh(p)
      zldis(p) = forc_hgt_u_pft(p) - displa(p)

      ! Check to see if the forcing height is below the canopy height
      if (zldis(p) < 0._r8) then
         found = .true.
         index = p
      end if

   end do

   if (found) then
      write(iulog,*)'Error: Forcing height is below canopy height for pft index ',index
      call endrun()
   end if

   do f = 1, fn
      p = filterp(f)
      c = pcolumn(p)

      ! Initialize Monin-Obukhov length and wind speed

      call MoninObukIni(ur(p), thv(c), dthv(p), zldis(p), z0mv(p), um(p), obu(p))

   end do

   ! Set counter for leaf temperature iteration (itlef)

   itlef = 0    
   fnorig = fn
   fporig(1:fn) = filterp(1:fn)

   ! Make copies so that array sections are not passed in function calls to friction velocity
   
   do f = 1, fn
      p = filterp(f)
      displa_loc(p) = displa(p)
      z0mv_loc(p) = z0mv(p)
      z0hv_loc(p) = z0hv(p)
      z0qv_loc(p) = z0qv(p)
   end do

   ! Begin stability iteration

   ITERATION : do while (itlef <= itmax .and. fn > 0)

      ! Determine friction velocity, and potential temperature and humidity
      ! profiles of the surface boundary layer

      call FrictionVelocity (lbp, ubp, fn, filterp, &
                             displa_loc, z0mv_loc, z0hv_loc, z0qv_loc, &
                             obu, itlef+1, ur, um, ustar, &
                             temp1, temp2, temp12m, temp22m, fm)

      do f = 1, fn
         p = filterp(f)
         c = pcolumn(p)
         g = pgridcell(p)

         tlbef(p) = t_veg(p)
         del2(p) = del(p)

         ! Determine aerodynamic resistances

         ram1(p)  = 1._r8/(ustar(p)*ustar(p)/um(p))
         rah(p,1) = 1._r8/(temp1(p)*ustar(p))
         raw(p,1) = 1._r8/(temp2(p)*ustar(p))

         ! Bulk boundary layer resistance of leaves

         uaf(p) = um(p)*sqrt( 1._r8/(ram1(p)*um(p)) )
         cf  = 0.01_r8/(sqrt(uaf(p))*sqrt(dleaf(ivt(p))))
         rb(p)  = 1._r8/(cf*uaf(p))
         rb1(p) = rb(p)
  
         ! Parameterization for variation of csoilc with canopy density from
         ! X. Zeng, University of Arizona

         w = exp(-(elai(p)+esai(p)))

         ! changed by K.Sakaguchi from here
         ! transfer coefficient over bare soil is changed to a local variable
         ! just for readability of the code (from line 680)
         csoilb = (vkc/(0.13_r8*(z0mg(c)*uaf(p)/1.5e-5_r8)**0.45_r8))

         !compute the stability parameter for ricsoilc  ("S" in Sakaguchi&Zeng,2008)

         ri = ( grav*htop(p) * (taf(p) - t_grnd(c)) ) / (taf(p) * uaf(p) **2.00_r8)

         !! modify csoilc value (0.004) if the under-canopy is in stable condition

         if ( (taf(p) - t_grnd(c) ) > 0._r8) then
               ! decrease the value of csoilc by dividing it with (1+gamma*min(S, 10.0))
               ! ria ("gmanna" in Sakaguchi&Zeng, 2008) is a constant (=0.5)
               ricsoilc = csoilc / (1.00_r8 + ria*min( ri, 10.0_r8) )
               csoilcn = csoilb*w + ricsoilc*(1._r8-w)
         else
              csoilcn = csoilb*w + csoilc*(1._r8-w)
         end if

         !! Sakaguchi changes for stability formulation ends here

         rah(p,2) = 1._r8/(csoilcn*uaf(p))
         raw(p,2) = rah(p,2)

         ! Stomatal resistances for sunlit and shaded fractions of canopy.
         ! Done each iteration to account for differences in eah, tv.

         svpts(p) = el(p)                         ! pa
         eah(p) = forc_pbot(g) * qaf(p) / 0.622_r8   ! pa
      end do

      ! 4/25/05, PET: Now calling the sun/shade version of Stomata by default
      call Stomata (fn, filterp, lbp, ubp, svpts, eah, o2, co2, rb, dayl_factor, phase='sun')
      call Stomata (fn, filterp, lbp, ubp, svpts, eah, o2, co2, rb, dayl_factor, phase='sha')

      do f = 1, fn
         p = filterp(f)
         c = pcolumn(p)
         g = pgridcell(p)

         ! Sensible heat conductance for air, leaf and ground
         ! Moved the original subroutine in-line...

         wta    = 1._r8/rah(p,1)             ! air
         wtl    = (elai(p)+esai(p))/rb(p)    ! leaf
         wtg(p) = 1._r8/rah(p,2)             ! ground
         wtshi  = 1._r8/(wta+wtl+wtg(p))

         wtl0(p) = wtl*wtshi         ! leaf
         wtg0    = wtg(p)*wtshi      ! ground
         wta0(p) = wta*wtshi         ! air

         wtga    = wta0(p)+wtg0      ! ground + air
         wtal(p) = wta0(p)+wtl0(p)   ! air + leaf

         ! Fraction of potential evaporation from leaf

         if (fdry(p) .gt. 0._r8) then
            rppdry  = fdry(p)*rb(p)*(laisun(p)/(rb(p)+rssun(p)) + &
                                     laisha(p)/(rb(p)+rssha(p)))/elai(p)
         else
            rppdry = 0._r8
         end if
         efpot = forc_rho(g)*wtl*(qsatl(p)-qaf(p))

         if (efpot > 0._r8) then
            if (btran(p) > btran0) then
               qflx_tran_veg(p) = efpot*rppdry
               rpp = rppdry + fwet(p)
            else
               !No transpiration if btran below 1.e-10
               rpp = fwet(p)
               qflx_tran_veg(p) = 0._r8
            end if
            !Check total evapotranspiration from leaves
            rpp = min(rpp, (qflx_tran_veg(p)+h2ocan(p)/dtime)/efpot)
         else
            !No transpiration if potential evaporation less than zero
            rpp = 1._r8
            qflx_tran_veg(p) = 0._r8
         end if

         ! Update conductances for changes in rpp
         ! Latent heat conductances for ground and leaf.
         ! Air has same conductance for both sensible and latent heat.
         ! Moved the original subroutine in-line...

         wtaq    = frac_veg_nosno(p)/raw(p,1)                        ! air
         wtlq    = frac_veg_nosno(p)*(elai(p)+esai(p))/rb(p) * rpp   ! leaf

         !Litter layer resistance. Added by K.Sakaguchi
         snowdp_c = z_dl ! critical depth for 100% litter burial by snow (=litter thickness)
         fsno_dl = snowdp(c)/snowdp_c    ! effective snow cover for (dry)plant litter
         elai_dl = lai_dl*(1._r8 - min(fsno_dl,1._r8)) ! exposed (dry)litter area index
         rdl = ( 1._r8 - exp(-elai_dl) ) / ( 0.004_r8*uaf(p)) ! dry litter layer resistance

         ! add litter resistance and Lee and Pielke 1992 beta
         if (delq(p) .lt. 0._r8) then  !dew. Do not apply beta for negative flux (follow old rsoil)
            wtgq(p) = frac_veg_nosno(p)/(raw(p,2)+rdl)
         else
            wtgq(p) = soilbeta(c)*frac_veg_nosno(p)/(raw(p,2)+rdl)
         end if

         wtsqi   = 1._r8/(wtaq+wtlq+wtgq(p))

         wtgq0    = wtgq(p)*wtsqi      ! ground
         wtlq0(p) = wtlq*wtsqi         ! leaf
         wtaq0(p) = wtaq*wtsqi         ! air

         wtgaq    = wtaq0(p)+wtgq0     ! air + ground
         wtalq(p) = wtaq0(p)+wtlq0(p)  ! air + leaf

         dc1 = forc_rho(g)*cpair*wtl
         dc2 = hvap*forc_rho(g)*wtlq

         efsh   = dc1*(wtga*t_veg(p)-wtg0*t_grnd(c)-wta0(p)*thm(p))
         efe(p) = dc2*(wtgaq*qsatl(p)-wtgq0*qg(c)-wtaq0(p)*forc_q(g))

         ! Evaporation flux from foliage

         erre = 0._r8
         if (efe(p)*efeb(p) < 0._r8) then
            efeold = efe(p)
            efe(p)  = 0.1_r8*efeold
            erre = efe(p) - efeold
         end if
         dt_veg(p) = (sabv(p) + air(p) + bir(p)*t_veg(p)**4 + &
              cir(p)*t_grnd(c)**4 - efsh - efe(p)) / &
              (- 4._r8*bir(p)*t_veg(p)**3 +dc1*wtga +dc2*wtgaq*qsatldT(p))
         t_veg(p) = tlbef(p) + dt_veg(p)
         dels = dt_veg(p)
         del(p)  = abs(dels)
         err(p) = 0._r8
         if (del(p) > delmax) then
            dt_veg(p) = delmax*dels/del(p)
            t_veg(p) = tlbef(p) + dt_veg(p)
            err(p) = sabv(p) + air(p) + bir(p)*tlbef(p)**3*(tlbef(p) + &
                 4._r8*dt_veg(p)) + cir(p)*t_grnd(c)**4 - &
                 (efsh + dc1*wtga*dt_veg(p)) - (efe(p) + &
                 dc2*wtgaq*qsatldT(p)*dt_veg(p))
         end if

         ! Fluxes from leaves to canopy space
         ! "efe" was limited as its sign changes frequently.  This limit may
         ! result in an imbalance in "hvap*qflx_evap_veg" and
         ! "efe + dc2*wtgaq*qsatdt_veg"

         efpot = forc_rho(g)*wtl*(wtgaq*(qsatl(p)+qsatldT(p)*dt_veg(p)) &
            -wtgq0*qg(c)-wtaq0(p)*forc_q(g))
         qflx_evap_veg(p) = rpp*efpot
         
         ! Calculation of evaporative potentials (efpot) and
         ! interception losses; flux in kg m**-2 s-1.  ecidif
         ! holds the excess energy if all intercepted water is evaporated
         ! during the timestep.  This energy is later added to the
         ! sensible heat flux.

         ecidif = 0._r8
         if (efpot > 0._r8 .and. btran(p) > btran0) then
            qflx_tran_veg(p) = efpot*rppdry
         else
            qflx_tran_veg(p) = 0._r8
         end if
         ecidif = max(0._r8, qflx_evap_veg(p)-qflx_tran_veg(p)-h2ocan(p)/dtime)
         qflx_evap_veg(p) = min(qflx_evap_veg(p),qflx_tran_veg(p)+h2ocan(p)/dtime)

         ! The energy loss due to above two limits is added to
         ! the sensible heat flux.

         eflx_sh_veg(p) = efsh + dc1*wtga*dt_veg(p) + err(p) + erre + hvap*ecidif

         ! Re-calculate saturated vapor pressure, specific humidity, and their
         ! derivatives at the leaf surface

         call QSat(t_veg(p), forc_pbot(g), el(p), deldT, qsatl(p), qsatldT(p))

         ! Update vegetation/ground surface temperature, canopy air
         ! temperature, canopy vapor pressure, aerodynamic temperature, and
         ! Monin-Obukhov stability parameter for next iteration.

         taf(p) = wtg0*t_grnd(c) + wta0(p)*thm(p) + wtl0(p)*t_veg(p)
         qaf(p) = wtlq0(p)*qsatl(p) + wtgq0*qg(c) + forc_q(g)*wtaq0(p)

         ! Update Monin-Obukhov length and wind speed including the
         ! stability effect

         dth(p) = thm(p)-taf(p)
         dqh(p) = forc_q(g)-qaf(p)
         delq(p) = wtalq(p)*qg(c)-wtlq0(p)*qsatl(p)-wtaq0(p)*forc_q(g)

         tstar = temp1(p)*dth(p)
         qstar = temp2(p)*dqh(p)

         thvstar = tstar*(1._r8+0.61_r8*forc_q(g)) + 0.61_r8*forc_th(g)*qstar
         zeta = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c))

         if (zeta >= 0._r8) then     !stable
            zeta = min(2._r8,max(zeta,0.01_r8))
            um(p) = max(ur(p),0.1_r8)
         else                     !unstable
            zeta = max(-100._r8,min(zeta,-0.01_r8))
            wc = beta*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8
            um(p) = sqrt(ur(p)*ur(p)+wc*wc)
         end if
         obu(p) = zldis(p)/zeta

         if (obuold(p)*obu(p) < 0._r8) nmozsgn(p) = nmozsgn(p)+1
         if (nmozsgn(p) >= 4) obu(p) = zldis(p)/(-0.01_r8)
         obuold(p) = obu(p)

      end do   ! end of filtered pft loop

      ! Test for convergence

      itlef = itlef+1
      if (itlef > itmin) then
         do f = 1, fn
            p = filterp(f)
            dele(p) = abs(efe(p)-efeb(p))
            efeb(p) = efe(p)
            det(p)  = max(del(p),del2(p))
         end do
         fnold = fn
         fn = 0
         do f = 1, fnold
            p = filterp(f)
            if (.not. (det(p) < dtmin .and. dele(p) < dlemin)) then
               fn = fn + 1
               filterp(fn) = p
            end if
         end do
      end if

   end do ITERATION     ! End stability iteration

   fn = fnorig
   filterp(1:fn) = fporig(1:fn)

   do f = 1, fn
      p = filterp(f)
      c = pcolumn(p)
      g = pgridcell(p)

      ! Energy balance check in canopy

      err(p) = sabv(p) + air(p) + bir(p)*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) &
         + cir(p)*t_grnd(c)**4 - eflx_sh_veg(p) - hvap*qflx_evap_veg(p)

      ! Fluxes from ground to canopy space

      delt    = wtal(p)*t_grnd(c)-wtl0(p)*t_veg(p)-wta0(p)*thm(p)
      taux(p) = -forc_rho(g)*forc_u(g)/ram1(p)
      tauy(p) = -forc_rho(g)*forc_v(g)/ram1(p)
      eflx_sh_grnd(p) = cpair*forc_rho(g)*wtg(p)*delt
      qflx_evap_soi(p) = forc_rho(g)*wtgq(p)*delq(p)

      ! 2 m height air temperature

      t_ref2m(p) = thm(p) + temp1(p)*dth(p)*(1._r8/temp12m(p) - 1._r8/temp1(p))
      t_ref2m_r(p) = t_ref2m(p)

      ! 2 m height specific humidity

      q_ref2m(p) = forc_q(g) + temp2(p)*dqh(p)*(1._r8/temp22m(p) - 1._r8/temp2(p))

      ! 2 m height relative humidity

      call QSat(t_ref2m(p), forc_pbot(g), e_ref2m, de2mdT, qsat_ref2m, dqsat2mdT)
      rh_ref2m(p) = min(100._r8, q_ref2m(p) / qsat_ref2m * 100._r8)
      rh_ref2m_r(p) = rh_ref2m(p)

      ! Downward longwave radiation below the canopy

      dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(g) + &
         emv(p)*emg(c)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))

      ! Upward longwave radiation above the canopy

      ulrad(p) = ((1._r8-emg(c))*(1._r8-emv(p))*(1._r8-emv(p))*forc_lwrad(g) &
         + emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb*tlbef(p)**3*(tlbef(p) + &
         4._r8*dt_veg(p)) + emg(c)*(1._r8-emv(p))*sb*t_grnd(c)**4)

      ! Derivative of soil energy flux with respect to soil temperature

      cgrnds(p) = cgrnds(p) + cpair*forc_rho(g)*wtg(p)*wtal(p)
      cgrndl(p) = cgrndl(p) + forc_rho(g)*wtgq(p)*wtalq(p)*dqgdT(c)
      cgrnd(p)  = cgrnds(p) + cgrndl(p)*htvp(c)

      ! Update dew accumulation (kg/m2)
      
      h2ocan(p) = max(0._r8,h2ocan(p)+(qflx_tran_veg(p)-qflx_evap_veg(p))*dtime)
      
      ! total photosynthesis

      fpsn(p) = psnsun(p)*laisun(p) + psnsha(p)*laisha(p)
      
#if (defined CN) && (defined C13)
      ! 4/14/05: PET
      ! Adding isotope code
      rc13_canair(p) = c13o2(p)/(co2(p)-c13o2(p))
      rc13_psnsun(p) = rc13_canair(p)/alphapsnsun(p)
      rc13_psnsha(p) = rc13_canair(p)/alphapsnsha(p)
      c13_psnsun(p) = psnsun(p) * (rc13_psnsun(p)/(1._r8+rc13_psnsun(p)))
      c13_psnsha(p) = psnsha(p) * (rc13_psnsha(p)/(1._r8+rc13_psnsha(p)))
      !write(iulog,*) p,ivt(p),btran(p),psnsun(p),psnsha(p),alphapsnsun(p),alphapsnsha(p)
#endif     
      
   end do

   ! Filter out pfts which have small energy balance errors; report others

   fnold = fn
   fn = 0
   do f = 1, fnold
      p = filterp(f)
      if (abs(err(p)) > 0.1_r8) then
         fn = fn + 1
         filterp(fn) = p
      end if
   end do

   do f = 1, fn
      p = filterp(f)
      write(iulog,*) 'energy balance in canopy ',p,', err=',err(p)
   end do

   end subroutine CanopyFluxes

!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: Stomata
!
! !INTERFACE:
   subroutine Stomata (fn, filterp, lbp, ubp, ei, ea, o2, co2, rb, dayl_factor, phase)
!
! !DESCRIPTION: 
! Leaf stomatal resistance and leaf photosynthesis. Modifications for CN code.

! !REVISION HISTORY:
! 22 January 2004: Created by Peter Thornton
! 4/14/05: Peter Thornton: Converted Ci from local variable to pps struct member
!    now returns cisun or cisha per pft as implicit output argument.
!    Also sets alphapsnsun and alphapsnsha.
! 4/25/05, Peter Thornton: Adopted as the default code for CLM, together with
!   modifications for sun/shade canopy.  Renamed from StomataCN to Stomata,
!   and eliminating the older Stomata subroutine
! 3/6/09: Peter Thornton; added dayl_factor control on Vcmax, from Bill Bauerle

! !USES:
     use shr_kind_mod , only : r8 => shr_kind_r8
     use shr_const_mod, only : SHR_CONST_TKFRZ, SHR_CONST_RGAS
     use clmtype
     use clm_atmlnd   , only : clm_a2l
     use spmdMod      , only: masterproc
     use pftvarcon    , only : nbrdlf_dcd_tmp_shrub
     use pftvarcon    , only : nsoybean, npcropmin
!
! !ARGUMENTS:
     implicit none
     integer , intent(in)    :: fn                 ! size of pft filter
     integer , intent(in)    :: filterp(fn)        ! pft filter
     integer , intent(in)    :: lbp, ubp           ! pft bounds
     real(r8), intent(in)    :: ei(lbp:ubp)        ! vapor pressure inside leaf (sat vapor press at tl) (pa)
     real(r8), intent(in)    :: ea(lbp:ubp)        ! vapor pressure of canopy air (pa)
     real(r8), intent(in)    :: o2(lbp:ubp)        ! atmospheric o2 concentration (pa)
     real(r8), intent(in)    :: co2(lbp:ubp)       ! atmospheric co2 concentration (pa)
     real(r8), intent(inout) :: rb(lbp:ubp)        ! boundary layer resistance (s/m)
     real(r8), intent(in)    :: dayl_factor(lbp:ubp) ! scalar (0-1) for daylength
     character(len=*), intent(in) :: phase         ! 'sun' or 'sha'
!
! !CALLED FROM:
! subroutine CanopyFluxes in this module
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in variables
! new ecophys variables (leafcn, flnr) added 1/26/04
!
     integer , pointer :: pcolumn(:)     ! pft's column index
     integer , pointer :: pgridcell(:)   ! pft's gridcell index
     integer , pointer :: ivt(:)         ! pft vegetation type
     real(r8), pointer :: qe25(:)        ! quantum efficiency at 25C (umol CO2 / umol photon)
     real(r8), pointer :: c3psn(:)       ! photosynthetic pathway: 0. = c4, 1. = c3
     real(r8), pointer :: mp(:)          ! slope of conductance-to-photosynthesis relationship
     real(r8), pointer :: tgcm(:)        ! air temperature at agcm reference height (kelvin)
     real(r8), pointer :: forc_pbot(:)   ! atmospheric pressure (Pa)
     real(r8), pointer :: tl(:)          ! leaf temperature (Kelvin)
     real(r8), pointer :: btran(:)       ! soil water transpiration factor (0 to 1)
     real(r8), pointer :: apar(:)        ! par absorbed per unit lai (w/m**2)
     real(r8), pointer :: leafcn(:)      ! leaf C:N (gC/gN)
     real(r8), pointer :: flnr(:)        ! fraction of leaf N in the Rubisco enzyme (gN Rubisco / gN leaf)
     real(r8), pointer :: sla(:)         ! specific leaf area, projected area basis (m^2/gC)
     real(r8), pointer :: fnitr(:)       ! foliage nitrogen limitation factor (-)
!
! local pointers to implicit inout variables
!
     real(r8), pointer :: rs(:)          ! leaf stomatal resistance (s/m)
     real(r8), pointer :: psn(:)         ! foliage photosynthesis (umol co2 /m**2/ s) [always +]
     real(r8), pointer :: ci(:)          ! intracellular leaf CO2 (Pa)
#if (defined C13)
     real(r8), pointer :: alphapsn(:)    ! 13C fractionation factor for PSN ()
#endif

!
! local pointers to implicit out variables
!
     real(r8), pointer :: lnc(:)         ! leaf N concentration per unit projected LAI (gN leaf/m^2)
     real(r8), pointer :: vcmx(:)        ! maximum rate of carboxylation (umol co2/m**2/s)
!
!
! !LOCAL VARIABLES:
!EOP
!
     real(r8), parameter :: mpe = 1.e-6_r8   ! prevents overflow error if division by zero
     integer , parameter :: niter = 3     ! number of iterations
     integer  :: f,p,c,g ! indices
     integer  :: iter    ! iteration index
     real(r8) :: ab      ! used in statement functions
     real(r8) :: bc      ! used in statement functions
     real(r8) :: f1      ! generic temperature response (statement function)
     real(r8) :: f2      ! generic temperature inhibition (statement function)
     real(r8) :: tc      ! leaf temperature (degree celsius)
     real(r8) :: cs      ! co2 concentration at leaf surface (pa)
     real(r8) :: kc      ! co2 michaelis-menten constant (pa)
     real(r8) :: ko      ! o2 michaelis-menten constant (pa)
     real(r8) :: atmp    ! intermediate calculations for rs
     real(r8) :: btmp    ! intermediate calculations for rs
     real(r8) :: ctmp    ! intermediate calculations for rs
     real(r8) :: q       ! intermediate calculations for rs
     real(r8) :: r1,r2   ! roots for rs
     real(r8) :: ppf     ! absorb photosynthetic photon flux (umol photons/m**2/s)
     real(r8) :: wc      ! rubisco limited photosynthesis (umol co2/m**2/s)
     real(r8) :: wj      ! light limited photosynthesis (umol co2/m**2/s)
     real(r8) :: we      ! export limited photosynthesis (umol co2/m**2/s)
     real(r8) :: cp      ! co2 compensation point (pa)
     real(r8) :: awc     ! intermediate calcuation for wc
     real(r8) :: j       ! electron transport (umol co2/m**2/s)
     real(r8) :: cea     ! constrain ea or else model blows up
     real(r8) :: cf      ! s m**2/umol -> s/m
     real(r8) :: rsmax0  ! maximum stomatal resistance [s/m]
     real(r8) :: kc25    ! co2 michaelis-menten constant at 25c (pa)
     real(r8) :: akc     ! q10 for kc25
     real(r8) :: ko25    ! o2 michaelis-menten constant at 25c (pa)
     real(r8) :: ako     ! q10 for ko25
     real(r8) :: bp      ! minimum leaf conductance (umol/m**2/s)
     ! additional variables for new treatment of Vcmax, Peter Thornton, 1/26/04
     real(r8) :: act25   ! (umol/mgRubisco/min) Rubisco activity at 25 C
     real(r8) :: act     ! (umol/mgRubisco/min) Rubisco activity
     real(r8) :: q10act  ! (DIM) Q_10 for Rubisco activity
     real(r8) :: fnr     ! (gRubisco/gN in Rubisco)
!------------------------------------------------------------------------------

     ! Set statement functions

     f1(ab,bc) = ab**((bc-25._r8)/10._r8)
     f2(ab) = 1._r8 + exp((-2.2e05_r8+710._r8*(ab+SHR_CONST_TKFRZ))/(SHR_CONST_RGAS*0.001_r8*(ab+SHR_CONST_TKFRZ)))

     ! Assign local pointers to derived type members (pft-level)

     pcolumn   => clm3%g%l%c%p%column
     pgridcell => clm3%g%l%c%p%gridcell
     ivt       => clm3%g%l%c%p%itype
     tl        => clm3%g%l%c%p%pes%t_veg
     btran     => clm3%g%l%c%p%pps%btran
     if (phase == 'sun') then
        apar   => clm3%g%l%c%p%pef%parsun
        rs     => clm3%g%l%c%p%pps%rssun
        psn    => clm3%g%l%c%p%pcf%psnsun
        ci     => clm3%g%l%c%p%pps%cisun
#if (defined C13)
        alphapsn  => clm3%g%l%c%p%pps%alphapsnsun
#endif
        sla    => clm3%g%l%c%p%pps%slasun
        lnc    => clm3%g%l%c%p%pps%lncsun   
        vcmx   => clm3%g%l%c%p%pps%vcmxsun   
     else if (phase == 'sha') then
        apar   => clm3%g%l%c%p%pef%parsha
        rs     => clm3%g%l%c%p%pps%rssha
        psn    => clm3%g%l%c%p%pcf%psnsha
        ci     => clm3%g%l%c%p%pps%cisha
        sla    => clm3%g%l%c%p%pps%slasha   
#if (defined C13)
        alphapsn  => clm3%g%l%c%p%pps%alphapsnsha
#endif
        lnc    => clm3%g%l%c%p%pps%lncsha   
        vcmx   => clm3%g%l%c%p%pps%vcmxsha
     end if

     ! Assign local pointers to derived type members (gridcell-level)

     forc_pbot => clm_a2l%forc_pbot

     ! Assign local pointers to derived type members (column-level)

     tgcm        => clm3%g%l%c%p%pes%thm

     ! Assign local pointers to pft constants
     ! new ecophys constants added 1/26/04

     qe25      => pftcon%qe25
     c3psn     => pftcon%c3psn
     mp        => pftcon%mp
     leafcn    => pftcon%leafcn
     flnr      => pftcon%flnr
     fnitr     => pftcon%fnitr

     ! Set constant values

     kc25  = 30._r8
     akc   = 2.1_r8
     ko25  = 30000._r8
     ako   = 1.2_r8
     bp    = 2000._r8

     ! New constants for CN code, added 1/26/04

     act25 = 3.6_r8
     q10act = 2.4_r8
     fnr = 7.16_r8
     
     ! Convert rubisco activity units from umol/mgRubisco/min -> umol/gRubisco/s

     act25 = act25 * 1000.0_r8 / 60.0_r8

     do f = 1, fn
        p = filterp(f)
        c = pcolumn(p)
        g = pgridcell(p)

        ! Initialize rs=rsmax and psn=0 because calculations are performed only
        ! when apar > 0, in which case rs <= rsmax and psn >= 0
        ! Set constants

        rsmax0 = 2.e4_r8
        cf = forc_pbot(g)/(SHR_CONST_RGAS*0.001_r8*tgcm(p))*1.e06_r8
        if (apar(p) <= 0._r8) then          ! night time
           rs(p) = min(rsmax0, 1._r8/bp * cf)
           psn(p) = 0._r8
           lnc(p) = 0._r8
           vcmx(p) = 0._r8
#if (defined C13)
           alphapsn(p) = 1._r8
#endif
        else                             ! day time
           tc = tl(p) - SHR_CONST_TKFRZ
           ppf = 4.6_r8 * apar(p)                  
           j = ppf * qe25(ivt(p))
           kc = kc25 * f1(akc,tc)       
           ko = ko25 * f1(ako,tc)
           awc = kc * (1._r8+o2(p)/ko)
           cp = 0.5_r8*kc/ko*o2(p)*0.21_r8

           ! Modification for shrubs proposed by X.D.Z
           ! Why does he prefer this line here instead of in subr.
           ! CanopyFluxes? (slevis)
           ! Equivalent modification for soy following AgroIBIS
#if (defined CNDV)
           if (ivt(p) == nbrdlf_dcd_tmp_shrub) btran(p) = min(1._r8, btran(p) * 3.33_r8)
#endif
           if (ivt(p) == nsoybean) btran(p) = min(1._r8, btran(p) * 1.25_r8)
           
           ! new calculations for vcmax, 1/26/04
           lnc(p) = 1._r8 / (sla(p) * leafcn(ivt(p)))
		   act = act25 * f1(q10act,tc)
#if (defined CN)
           if ( ivt(p) < npcropmin )then
              vcmx(p) = lnc(p) * flnr(ivt(p)) * fnr * act / f2(tc) * btran(p) * &
                        dayl_factor(p)
           else
              vcmx(p) = 101._r8 * f1(q10act,tc) / f2(tc) * btran(p) * dayl_factor(p)
           end if
#else
           vcmx(p) = lnc(p) * flnr(ivt(p)) * fnr * act / f2(tc) * btran(p) * &
                     dayl_factor(p) * fnitr(ivt(p))
#endif
           
           ! First guess ci

           ci(p) = 0.7_r8*co2(p)*c3psn(ivt(p)) + 0.4_r8*co2(p)*(1._r8-c3psn(ivt(p)))

           ! rb: s/m -> s m**2 / umol

           rb(p) = rb(p)/cf 

           ! Constrain ea

           cea = max(0.25_r8*ei(p)*c3psn(ivt(p))+0.40_r8*ei(p)*(1._r8-c3psn(ivt(p))), min(ea(p),ei(p)) ) 

           ! ci iteration for 'actual' photosynthesis

           do iter = 1, niter
              wj = max(ci(p)-cp,0._r8)*j/(ci(p)+2._r8*cp)*c3psn(ivt(p)) + j*(1._r8-c3psn(ivt(p)))
              wc = max(ci(p)-cp,0._r8)*vcmx(p)/(ci(p)+awc)*c3psn(ivt(p)) + vcmx(p)*(1._r8-c3psn(ivt(p)))
              we = 0.5_r8*vcmx(p)*c3psn(ivt(p)) + 4000._r8*vcmx(p)*ci(p)/forc_pbot(g)*(1._r8-c3psn(ivt(p))) 
              psn(p) = min(wj,wc,we) 
              cs = max( co2(p)-1.37_r8*rb(p)*forc_pbot(g)*psn(p), mpe )
              atmp = mp(ivt(p))*psn(p)*forc_pbot(g)*cea / (cs*ei(p)) + bp
              btmp = ( mp(ivt(p))*psn(p)*forc_pbot(g)/cs + bp ) * rb(p) - 1._r8
              ctmp = -rb(p)
              if (btmp >= 0._r8) then
                 q = -0.5_r8*( btmp + sqrt(btmp*btmp-4._r8*atmp*ctmp) )
              else
                 q = -0.5_r8*( btmp - sqrt(btmp*btmp-4._r8*atmp*ctmp) )
              end if
              r1 = q/atmp
              r2 = ctmp/q
              rs(p) = max(r1,r2)
              ci(p) = max( cs-psn(p)*forc_pbot(g)*1.65_r8*rs(p), 0._r8 )
           end do

           ! rs, rb:  s m**2 / umol -> s/m 

           rs(p) = min(rsmax0, rs(p)*cf)
           rb(p) = rb(p) * cf 
           
#if (defined C13)
           ! 4/14/05: PET
           ! Adding isotope code
           alphapsn(p) = 1._r8 + (((c3psn(ivt(p)) * (4.4_r8 + (22.6_r8*(ci(p)/co2(p))))) + &
                         ((1._r8 - c3psn(ivt(p))) * 4.4_r8))/1000._r8)
           !alphapsn(p) = 1._r8
           !write(iulog,*) 'in StomataCN ',p,ivt(p),c3psn(ivt(p)),ci(p),co2(p),alphapsn(p)
#endif
           
        end if

     end do

  end subroutine Stomata

end module CanopyFluxesMod
